function x = f9(blk)
%F9 extracts 9 features from blk
[M,N] = size(blk);
if mod(M, 2) ~= 0 || mod(N, 2) ~= 0
    error('M and N must be multiplies of 2.');
end

% calculate F
f = zeros(1, 9);
x = f;
M2 = M / 2;
N2 = N / 2;
f(1) = sum(blk(:)) / M / N;
f(2) = sum(sum(blk(1:M2, 1:N2)));
f(3) = sum(sum(blk(1:M2, N2+1:N)));
f(4) = sum(sum(blk(M2+1:M, 1:N2)));
f(5) = sum(sum(blk(M2+1:M, N2+1:N)));
f(6:9) = f(2:5) - f(1);
f(2:5) = f(2:5) / (4 * f(1) + eps);

m1 = max(f(6:9));
m2 = min(f(6:9));

x(1) = floor(f(1));
x(2:5) = floor(255 * f(2:5));
x(6:9) = floor(255 * (f(6:9) - m2) / (m1 - m2 + eps) );
end

function out = radixsort(v, M)
    if nargin < 2
        M = 255;
    end
    p = numel(v);
    h = histc(v,[0:M]);
    out = zeros(1, p);
    count = 0;
    for i = 0:M
        if h(i) > 0
           out(count+1:count+h(i)) = i;
           count = count + h(i);
        end
    end
end
